Construction of localized wave functions for a disordered optical lattice and analysis 

of the resulting Hubbard model parameters 
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We propose a method to construct localized single particle wave functions using imaginary time 
projection and thereby determine lattice Hamiltonian parameters. We apply the method to a specific 
disordered potential generated by an optical lattice experiment and calculate for each instance of 
disorder, the equivalent lattice model parameters. The probability distributions of the Hubbard 
parameters are then determined. Tests of localization and eigen-energy convergence are examined. 

PACS numbers: Valid PACS appear here 
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I. INTRODUCTION 

Understanding the properties of disordered materi- 
als has a fundamental significance in condensed matter 
physics. Various kinds of disorder exist in real materials, 
but their disorder is difficult to characterize and control 
experimentally. Recently developed optical lattice tech- 
niques [l| have enabled the construction of a nearly per- 
fectly controlled disordered potential and the measure- 
ment of properties of strongly correlated atoms in that 
potential provides an opportunity to compare quantita- 
tively experimental results with parameter-free theoreti- 
cal calculations. Bosonic atoms are particularly interest- 
ing because then Monte Carlo simulation is efficient. 

In this paper, we consider the problem of mapping a 
disordered single body potential to a lattice model. By 
removing the high energy states associated with the con- 
tinuum, the Monte Carlo simulation becomes more effi- 
cient. Particularly, efficient algorithms have been devel- 
oped [2] 3] for lattice models. Consider the continuum 
Hamiltonian of N atoms with mass m moving in the ex- 
ternal potential U(r) and interacting with the pairwise 
potential energy V(r a — rp) 
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where indices a, (3 label the atoms. On the other hand, 
the quantum mechanics of particles moving in a lattice is 
conveniently described in a basis of localized wave func- 
tions, such as the Wannier functions associated with a pe- 
riodic potential. Using these localized functions, we can 
define an effective lattice Hubbard Hamiltonian. Written 
in second quantized notation it has the form: 



h = - tija\aj + CjUj + ^ Ujrii{rii — 1) 

- E h a \ a i + Uijniiij + (2) 



(ij) 



where i labels the single particle states (lattice sites), (ij) 
denotes a nearest neighbor pair and {ij} a next-nearest 
neighbor pair, and tij are hopping coefficients, q is 
the on-site energy, m is the on-site interaction and «y 
is the nearest neighbor off-site interaction. (Terms such 
as next-nearest neighbor hopping and offsite interaction 
are often neglected.) Note that Hn refers to the TV-body 
Hamiltonian in continuous space and h to its equivalent 
on a lattice. 

In a periodic potential, Wannier functions of a given 
band are related to the Bloch functions ip n k of the same 
band n by the unitary transformation 

W ni (r)=w n (r-TL i ) = J=^„ k (r) e - ik ' R '. (3) 
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w n i is localized around the lattice site Ri However, in 
the absence of periodicity, the concept of Wannier func- 
tions needs to be generalized. Two main types of gener- 
alizations exist in literature. The perturbative approach 
5] assumes the existence of the band structure and thus 
is applicable to nearly periodic potentials. The varia- 
tional approach Q Q [8[ emphasizes the minimization of 
the spatial spread with respect to unitary transforma- 
tions of a starting basis set, for example the Wannier 
functions of a periodic potential. 

In order to be useful, we would like the generalized 
Wannier functions to have the following properties. First, 
localization is required by the physical picture of parti- 
cles hopping in the lattice. Second, a correct descrip- 
tion of the low energy density of states is necessary to 
capture the low temperature physics. Third, for conve- 
nience, the orthogonality of the basis set is required to 
use commutation relations of creation and annihilation 
operators in the second quantized Hamiltonian. Finally, 
we would like the lattice Hamiltonian to be free of the 
sign problem so that quantum Monte Carlo calculations 
are efficient. This requires the off-diagonal elements to 
be non-positive (i.e. > 0). Note that the original 
Hamiltonian 1-Ljq has this property. 

In section II, we propose a method of constructing lo- 
calized single particle basis functions based on imaginary 
time evolution of localized basis functions: Wi(0) where 
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labels the site. 
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where H.\ denotes the one particle continuum Hamilto- 
nian. This has the effect of suppressing the high energy 
components but also spreading out the basis states. In 
section III, results of the method are presented for the 
specific disorder probed experimentally. 



respect to r and multiply on the right and left by S 1 ^ 2 , 
we find an expression for h in terms of 5*: 



z_idS a_i 
-S 2 —S 2 
dr 



' e (*-*)M^j e (*-*)*dA. (8) 



If we assume that h becomes r-independent as r — > oo, 
we can neglect the second term on the right hand side 
and find: 



II. METHOD 

In constructing a lattice model, our goal is to coarse 
grain the description of the continuum system, so that 
instead of recording the precise position of an atom, we 
only record which lattice site it occupies. We match up 
the lattice and continuum models using the density ma- 
trix; we require that the low temperature density matrix 
of the lattice model to be identical to the reduced den- 
sity matrix of the continuum system when high energy 
degrees of freedoms are traced out. Use of the density 
matrix is motivated by the fact that the linear response 
of a system to an external perturbation, either an exter- 
nal field, or particle insertion, or a coupling to another 
subsystem is determined by its one-body density matrix 
@[T3|. If we match the density matrices, the lattice sys- 
tem is guaranteed to have not only the same density dis- 
tribution n(r) and hopping properties, such as diffusion, 
but also the same response to external perturbations as 
the continuum system. 

The unnormalized single particle density matrix in the 
continuum system is defined by: 



,r';r) 



-tHi 
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Let Wj(r;0) be a localized basis which assigns atoms to 
lattice sites, e.g. Wi(r;Q) = 1 if |r — Ri| is minimized 
with respect to i, i.e. in the i th Wigner-Seitz cell. Then 
a course-grained density matrix is defined as 



Sy(r) = ^(0)|e-^K-(0) 



drdr'w*(r; 0)p(r, r'; r)wj (r'; 0) (6) 



Note that if Wi(r;0)'s are chosen to be everywhere pos- 
itive, all elements of the lattice density matrix Sij are 
also positive and can be used in a lattice QMC calcula- 
tion directly. We now want to construct a single-particle 
Hamiltonian which when solved gives Sij(r) for large r, 
or in matrix notation to determine h such that 



S(r) 



-rh 



(7) 



[l6| Formally, the solution h = — r _1 logS 1 (t) may have 
some r dependance and not necessarily have the other 
properties mentioned above. Differentiating Eq. ([7]) with 



dr 



(9) 



Consider the eigenfunction expansion of the continuum 
density matrix: 



p(r,r';r)=E^( r )^( r ') e " Ti 



(10) 



where E a and 4> a are the 1-particle eigenvalues and eigen- 
functions of the continuum Hamiltonian. For a suffi- 
ciently large r, and for a system with a gap, only states 
below the gap will survive. If there are TV such states, 
it is clear that we will capture the density of states with 
exactly N basis functions Wi. Now let us define the or- 
thogonalized basis by splitting up the density operator 

exp (-tWi) = exp {-\rU^j exp (-\tH^J (11) 

and having it act partially to the left and right in Eq. ©. 
Combining Eq. ^ and Eq.® we obtain the expression 
for the model Hamiltonian: 

hij = E S ^(r) (w k {r 12)^x^/2)) S^{r) 



= («} i (T/2)|7l 1 |tS J (r/2) 



(12) 



where Wi(r) 



Wi(0) are the non-orthonormalized 



basis functions at time r and w(t/2) = S 1 / 2 (t)w(t/2) 
are the orthonormalized basis functions because 



S^t) = H(t/2)K(t/2)) 



(13) 



is the overlap matrix. This is known as Lowdin orthogo- 
nalization [ll|[T7j. 

The imaginary time evolution is equivalent to a dif- 
fusion process with sinks or sources determined by the 
potential U(r) . Without a potential present, an initially 
localized distribution will spread out as \fr as a func- 
tion of imaginary time. When the wavepacket (or basis 



function) \wi{t)) 



-tHi 



\wi(0)} encounters the regions 



of high potential energy separating the lattice sites, its 
diffusion will stop, until it tunnels through to the next 
site. If the assumption of temperature-independence 



(14) 
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is correct, according to Eq. (fT2]) . the orthogonalized basis 
§-1/2 ^2r) \wi(r)) converges at large r. Instead of taking 
the logarithm of the reduced density matrix Eq. |7|). we 
choose to construct the lattice Hamiltonian from Eq. © 
for two reasons. Firstly, numerical tests show that Eq. © 
converges faster than — ilogS* as r increases. Secondly, 
the explicit construction of basis functions enables us to 
calculate the interaction terms in the second quantized 
many body Hamiltonian. Finally, use of Eq. © instead 
of Eq. ((7|) gives a spectrum as upper bounds to the con- 
tinuum spectrum. 

The choice of the initial basis function Wi (r; 0) is to 
some extent arbitrary, as long as it is non-negative and 
localized within a lattice cell, and there is one basis func- 
tion for each lattice site. We choose to set Wi(r; 0) = 1 
inside a cube of side a centered on R,. 



second FFT, Eq. (fT7|) . we add a buffer layer outside the 
cube with thickness chosen so that the localization re- 
gion of the evolved function over one imaginary time step 
does not exceed the cube in which FFT is performed; the 
thickness is proportional to •v/r/m. We periodically ex- 
amine the evolved basis set, to determine if the cube can 
be made smaller. A common normalization factor is re- 
quired for all basis functions every several steps to avoid 
numerical overflow or underflow. 

Eq. © demands orthogonalization of the basis set. 
Lowdin orthogonalization preserves, as much as possi- 
ble, the localization and symmetry of the original non- 
orthogonal basis states. In terms of the overlap matrix, 
we construct a set of orthogonalized states 



= £(*- 1/2 )d 



(19) 



A. Algorithm for imaginary time projection and 
orthogonalization 

To apply the imaginary time evolution to the con- 
struction of localized wave functions and thereby to ex- 
tract microscopic parameters of the corresponding lattice 
model, we then start with N initial trial wave functions, 
each of which is well localized in one lattice cell. Each 
wave function is independently evolved over a sufficiently 
long imaginary time. The set of N wave functions are 
then transformed into an orthonormal basis. 

To perform the imaginary time evolution, consider the 
Trotter formula! 12] 
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(15) 



In a coordinate representation, a single step of imaginary 
time t can be written as: 

w{r,t + r) = J dV(r|e- r7il |rV(r',t) 



This integral is a convolution, and can be efficiently eval- 
uated by Fast Fourier Transform 



w(r, t + r) = FFT 



e 2 " 



7k 



(17) 



where /k is defined as an inverse- Fourier transform 



/k = FFT- 



~w(r, t) 



(18) 



We can also take advantage of the localization of w(r): 
it is vanishingly small away from its initial site, so that 
we only store its values in a cube enclosing the region 
in which the wave function is non-zero. When doing the 



No other set of orthonormal states generated from the 
space spanned by the original non-orthogonal set of states 
resemble the original set more closely, in the sense of least 
square, than do Lowdin set of states [3. Explicitly, 
Lowdin orthogonalization minimizes the expression 



N 



Wi\ 



(20) 



over all linear transformations T which transform the 
original non-orthonormal set of states \wi) into an or- 
thonormal set of states Im,-) 



{Wi,Wj) 



Twi,Tw, 



(21) 



For efficiency, this procedure can be done in an itera- 
tive fashion. Because the original non-orthogonal set of 
wave functions are localized, the overlap matrix S mn has 
the form of the identity matrix plus a small off-diagonal 
part 



■At 



(22) 



where the diagonal elements of A are zero and the off- 
diagonal elements |-A m „| <C 1. This enables us to per- 
form Lowdin orthogonalization iteratively by repeated 
application of the approximate inverse square root of the 
overlap matrix 



(s- 1/2 h 



(23) 



to the non-orthogonal basis set by updating the overlap 
matrix at each step 13]. Hence the basis set is iterated 



(24) 



until convergence is reached, \A\ ~ 0. The convergence 
of the overlap matrix to identity matrix is geometric. 
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The iterative scheme is efficient for large systems be- 
cause the basis sets are sparse. The computation time of 
this algorithm is linear in the number of lattice sites, i.e. 
the complexity is proportional to 



(# of steps) ■ M ■ n log i 



(25) 



where M is the number of lattice sites and n is the num- 
ber of pixels for each basis function. 



B. Hubbard parameters 

Once the orthogonalized basis set has been con- 
structed, the effective lattice Hamiltonian is obtained. 
For convenience we drop the r dependance. According 
to Eq. (|12p . the single particle Hubbard parameters are 
calculated as detailed in Eq. ([26]) - ((27| : the on-site ener- 
gies 



w*(r)Wi^(r)d 3 



and the hopping coefficients 



w* (r)Ti.iWj(r)d 3 r. 



(26) 



(27) 



The interaction term is computed from first-order pertur- 
bation theory in V. In the case of a contact interaction 
with the scattering length a s we find for u: 



4Tra s h 2 

Ui = 

m 

and the off-site interaction 



Wi(r)| z \ wj(r)\ 2 d 3 r 



(28) 



(29) 



The problem of a single particle moving in a periodic 
potential of the form cos x + cos y + cos z can be solved 
analytically (l4| . We compared the imaginary time pro- 
jected states with the results from exact diagonalization; 
this is shown in Fig. [1] and Fig. O We used a spatial grid 
with 8 3 pixels per lattice cell and an imaginary time step 



At = io- 4 £;; 



We find the error vanishes linearly as 



the time step goes to zero. 



C. Measures of localization and energy 

Several quantities can be used to characterize the local- 
ization and accuracy of the basis set. The spatial spread 
fl w = (r 2 ) w — (r)^ quantifies the localization of a wave 
function. 

The off-site integral Uij measures the spatial overlap 
between a pair of basis states. If it is small relative to tij 
and the Hubbard U, the approximation of keeping only 
the on-site interaction in the lattice model is appropriate. 
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FIG. 1: Hubbard U as a function of the potential depth of 
lattice field; the line is obtained from exact diagonalization, 
and the open circles are obtained using the method in this 
paper. 
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FIG. 2: Log-scale hopping coefficient as a function of the 
potential depth of the lattice field; the line is obtained from 
exact diagonalization and the open circles are obtained using 
the method in this paper. 



Its rms value over all nearest neighbor pairs measures the 
whole basis set. 

The convergence rate of the N x N matrix of the sin- 
gle particle lattice Hamiltonian is measured by the time 
derivatives of its N eigenvalues E^ ttice 's, 



r 



i 



\ " 



TV 4^ 



dT 



-E- 



(30) 



To determine the accuracy of the basis set we compare 
r "' ce 's with the lowest N eigenvalues E^ act of the 



original continuum Hamiltonian ri± estimated from 



E, 



(?) 



E- 



(0 



3 (t -> oo). 



(31) 
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The worst case error is defined as 



eia 



E [l) - E 



exact 



(32) 



We did the same estimate for the lattice Hamiltonian that 
has only nearest neighbor hopping terms: we denote this 



III. RESULTS FOR A DISORDERED LATTICE 

We now apply our method to the disordered lattice 
potential created in the White et. al. experiment [l|, 
87 Rb atoms are trapped in a background cubic lattice 
potential created by red lasers with wave vector k = —. 
The periodic potential is: 



Ul(t) 



-Sl x ^ cos 

1=1 



27rn 7 - 



(33) 



where n^, i = 1,2,3 are three mutually orthogonal unit 
vectors. A disordered speckle field £7o(r) is produced by 
a laser beam with phases randomized by a diffuser. The 
speckle field is everywhere positive, characterized by a 
speckle strength Sr>- 

U D (r) > 0, Vr 
(U D (r)) » 0.75S D 
(f/f,(r)) « S D 

and a spatial auto-correlation (Ud(^)Ud(^')) with corre- 
lation length ~ 1.29a, i.e. slightly larger than the lattice 
spacing. 

The total external potential is a superposition of 
the periodic lattice potential and the speckle potential 
U(r) = Ul(t) + Ud{y)- The single particle Hamiltonian 
in units of the recoil energy En = is: 



Hi = — 



V 2 



U(r). 



(34) 



We constructed our potential to match the experiment as 
closely as possible; see the reference p| for more details. 

To investigate the evolution of lattice Hamiltonian 
Eq. @, at every step of the imaginary time, the basis 
set is orthonormalized before constructing the Hamilto- 
nian matrix and calculating the energies E ] 



lattii 



Then 

the basis set is evolved using the previous basis set be- 
fore orthonormalization; this means each basis function 
is evolved independently. 

To illustrate the convergence of the matrix elements of 
the lattice Hamiltonian, the evolution diagram of an on- 
site energy on one particular site and a nearest neighbor 
hopping coefficient on one particular bond for Sl = 14 
and Sjy = 1 are shown in Fig. [3] We characterize the 
localization of the basis functions by the average nearest 
neighbor off-site integral «y, which measures the spatial 
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- Nearest Neighbor Hopping t 



0.0079 E_ 



0.01 LU 



Imaginary Time (E R 1 ) 

FIG. 3: Evolution diagram of an on-site energy(left scale) and 
a nearest neighbor hopping coefficient (right scale) in a lattice 
for Sl = 14 and Sd = 1. At large imaginary time r, these 
two matrix elements approach constant values. 



overlap between a pair of nearest neighbor basis func- 
tions. Figs.LHshows the evolution diagrams of the average 
on-site interaction m and the average off-site interaction 
Uij, which are also converging at large imaginary time. 
The limiting value of the off-site interaction is 4 ~ 5 
orders of magnitude smaller than that of the on-site in- 
teraction, which means that the basis functions are still 
localized at large imaginary time. Note that although 
the imaginary time projection operator e~ rWl spreads 
out the basis states, Lowdin orthogonalization operator 
5-1/2 res tores their localization form. 

To illustrate the effect of Lowdin orthogonalization on 
the localization property of the basis set, the evolution 
diagram for Sl = 14 and Sd = 1 is shown in Fig. [5] by 
including the off-site integral of the set before orthogo- 
nalization. It can be seen from the graph that Lowdin 
procedure helps to localize the basis functions w(t). 

The localization characterized by the spatial spread 
Fl w and drift D w is shown in Fig.LH The maximum value 
among all basis functions is chosen to measure the whole 
basis set. As shown in the graph, the values that these 
two quantities asymptotes to at large time are small com- 
pared to the lattice constant, which signifies that the ba- 
sis functions are localized. 

The convergence rate T of eigen-energies of the lattice 
Hamiltonian is shown in Figs. [7] It can be seen from 
the graph that the effective lattice Hamiltonian becomes 
temperature-independent at low temperature. It is also 
illuminating to look at the evolution diagram of the worst 
case error Eq. (|32|) . as shown in Figs. [8] where the exact 
eigen-energies are estimated by 
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^lattice V 
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We compared the worst case error in energy for the near- 
est neighbor model (f = 0) versus the full lattice model. 
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FIG. 4: Evolution diagram of the average on-site interaction 
«(left scale) and the average nearest neighbor off-site inter- 
action w(right scale) for Sl = 14 and Sd = 1. Note that 
u measures the localization of a pair of basis functions. The 
diagram shows that limiting value of u is 4 ~ 5 orders of mag- 
nitude smaller than that of u, which indicates that the basis 
functions at large imaginary time are still localized. 
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FIG. 6: Evolution diagram of the maximum spatial spread 
and drift (average deviation from the initial position) in units 
of the lattice constant for Sl = 14 and Sd = 1. The values 
that these two quantities approach at large imaginary time 
are small compared to the lattice constant a, which means 
that the localization the basis functions is preserved. 
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FIG. 7: Convergence rate F of eigen-energies for Sl = 14 and 
Sd — 1 shown in log-scale. 
1 - 



FIG. 5: The effect of Lowdin orthogonalization on the offsitc 
integral for Sl = 14 and Sd = 1.0. It can be seen from the 
graph that the imaginary time evolution spreads out the wave 
packets, but Lowdin orthogonalization restores the localiza- 
tion. 



The spatial overlap between basis functions remain neg- 
ligible at the early stage so that the nearest neighbor 
model has almost the same spectrum as the full lattice 
model; the error in energy is reduced as imaginary time 
evolves. A finite error persists in the eigen-energy of the 
nearest neighbor model because the next nearest neigh- 
bor hopping terms are neglected. Note that this error is 
less than 10~ 4 Er. 

To explain how it is possible to suppress the energy 
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FIG. 8: Imaginary time evolution of the worst case error in 
energy for Sl = 14 and Sd = 1, shown in log- scale. 
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Eigen-Energy (Er) 

FIG. 9: (Color online) Density of states for a single particle 
in a disordered lattice with Sl = 14 and Sd = 1, 2, 3. Energy 
intervals with low density of states still exist in the presence 
of disorder. 




On-site Energy (E R ) 

FIG. 10: (Color online) Probability distribution of the on- 
site energy for Sl = 14 and Sd ~ 1. The line is a fit to an 
exponential function. 



of the original localized basis set before causing signifi- 
cant derealization, it is useful to look at the single par- 
ticle density of states. In particular, we are interested in 
whether the gap between bands persists in the presence 
of disorder. Fig. [9] shows the density of states of a single 
particle in the disordered lattice. 15 samples each with 
a 5 3 lattice for each disorder strength were calculated. 
It can been seen from the plot that the lowest band is 
broadened and skewed by the disorder; there remains a 
minimum in the density of states between the first band 
and the second band (pseudo-gap). It is the existence of 
such a gap that enables us to project out the high en- 
ergy components in the initial set of trial states before 
delocalization sets in. 




Nearest Neighbor Hopping (E 



IV. STATISTICS OF HUBBARD PARAMETERS 

We now discuss the statistical properties of the calcu- 
lated Hubbard parameters. These are shown in Figs. ITD1 
-[Hfor S L = 14 and S D = 1. About 1000 samples of 6 3 
sites are used to construct the probability distributions 
of Hubbard parameters. 

Fig. [TU] shows the probability distribution of the on- 
site energy e$ for Sd — 1 and Sl = 14. The distribution 
is skewed with a steep onset at low energy and a tail at 
high energy. We fit the distribution to an exponential 
decay function 



P(e) ~ exp (-e/r) 



(35) 



for s > -10.5E R finding T w 0.97 E R for S D = 1 and 
Sl = 14. Note that the disorder potential is always pos- 
itive, so that the on site energy is greater than its value 
for the periodic lattice which is — 10.58.Er for this value 
of Si. 

Hopping coefficients tij characterize the mobility of the 
atoms. Recall that negative values of t will cause dif- 
ficulty in quantum Monte Carlo calculations. Fig. [11] 



FIG. 11: Probability distribution of the nearest neighbor hop- 
ping with Sl = 14 and Sd = 1. This is a predominantly 
positive asymmetric distribution. 



shows the probability distribution of nearest neighbor 
hopping coefficients. This distribution is asymmetrically 
centered around its value 8 x 10~ 3 E R for the periodic 
potential with a width 



St 



= 0.15 



(36) 



In 10 6 sampled bonds, only positive t^ were found. 
For Sl = 14 and for 0.05 < Sd < 1, St/ (t) ranges from 
10~ 2 to 10" 1 . 

Fig. [12] shows the probability distribution of next- 
nearest neighbor hopping coefficients. This distribu- 
tion is symmetrically centered around zero with a width 
W = 1.25 x 10 _5 -Er and about 2 orders of magnitude 
smaller than nearest neighbor hopping. Note that in the 
clean limit, the next nearest neighbor hopping coefficient 
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Next Nearest Neighbor Hopping (10~ 5 E R ) 

FIG. 12: Probability distribution of the next nearest neigh- 
bor hopping for Sl = 14 and Sd = 1. This distribution is 
symmetrically centered around zero. 
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On-site Interaction (E R ) 

FIG. 13: (Color online) Probability distribution of the on-site 
interaction, i.e. Hubbard U, for Sl = 14 and Sd ~ 1. The 
line is a fit to a Laplace function. 
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is exactly zero for the simple cubic lattice by symmetry. 
As shown in Fig. [8j setting t = changes the resulting 
single particle energies by a maximum of 2 x 10~ 5 Er. 

Fig- ESI shows the probability distribution of the Hub- 
bard U, which is characterized by its narrow peak roughly 
centered around the value of u in the periodic limit 
(0.36.%.) with a 1% relative width. We fit this distri- 
bution to Laplace function 

P(u) = — exp _ j (37) 

with u w 0.36% and A = 2 x 10" 3 %. For S L = 14 
and for 10~ 2 < Sd < 1, 5u/ (u) ranges from 10~ 4 to 
10~ 2 . Hence one can assume that the on-site interac- 
tion is roughly constant even in the presence of disor- 
der. Fig. [14] shows the probability distribution of near- 
est neighbor overlap u. We observe that the magnitude 
of off-site interactions is almost 4 orders of magnitude 
smaller than the on-site interaction; evidently negligible 
in the many-body interacting Hamiltonian. 

On-site energies are usually assumed to be almost un- 
corrected between different sites. We calculated the 
nearest neighbor covariance function. For Sl — 14 and 
Sd = 1, with (ij) nearest neighbor pairs: 

(HEiLzMpl „ o.49 (38) 
(a*) -(e) 2 

The £j's are correlated between nearest neighboring sites 
for this disordered potential. 

Fig. [T5] shows the correlation pattern between the on- 
site energy difference of nearest neighboring sites and the 
hopping coefficient. Fit to this joint distribution gives 
- *o ~ (N ~ tj\) S with S = 1.05. In White et. 
al.Q], the orientation of laser speckles does not coin- 
cide with the lattice axes; the laser speckle has a cylin- 
drical shape whose longitudinal direction points along 
ini + |n 2 + 4g n 3 . The insert of Fig. [15] displays the 
correlation pattern for bonds in ri3-direction if the longi- 
tudinal direction of the speckles is aligned with the 113- 
axis of the lattice. This illustrates the directional effect 
of laser speckles. 

The characteristics of the speckle field is reflected in 
other aspects of the parameters. According to the orien- 
tation of laser speckles with respect to the lattice axes, we 
should expect that the average hopping coefficient along 
ni and ri2 to be equal and the hopping along 113 to be 
different. As shown in Table fl] (t z ) differs from those of 
(tx,y) because of the cylindrical symmetry of the speckle. 
However, the difference is small because the correlation 
length of the speckle is only slightly larger than the lat- 
tice spacing, such that the nearest neighbor hopping is 
not sensitive to the anisotropy induced by the speckle. 



FIG. 14: Probability distribution of the nearest neighbor off- 
site interaction for Sl ~ 14 and Sd ~ L 



In Fig.rjUthe variation of the distribution widths of all 
4 Hubbard parameters versus speckle intensity is shown 
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FIG. 16: (Color online) The width of the probability distri- 
bution for e, tij, Ui and (it) for Sl = 14. 



FIG. 15: Correlation between the on-site energy difference 
and hopping coefficient between nearest neighbor sites for 
Sl = 14 and Sd = 1- The insert displays the correlation 
pattern for bonds in ri3-direction if the longitudinal direction 
of the speckles is aligned with the 113-axis of the lattice. The 
line in the insert is a fit to a power function. 



TABLE I: Directional effect of speckles for Sl = 14 



Sp{E R ) {t x )(lQ- A E R ) fa)(l(T J £fl) {t z )(W-- i E R ) 

0.050 8.00 ±2xlCr 4 8.00 ±2xlCT 4 8.00 ± 2 x lCT* 

0.100 8.02 ±4 xl0~ 4 8.02 ±4 xlO~ 4 8.01 + 3x10^ 

0.250 8.10 ^1x10^ 8.10±lxl0~ 3 8.07+1x10^ 

0.375 8.20±2xl0~ a 8.20 ± 2 x 10~ 3 8.16 + 1x10^" 

0.500 8.32 ±3xl0~ a 8.33 ± 3 x 10~ 3 8.26 + 2x10^" 

0.750 8.59 ±4xl0~ 3 8.60 ± 4 x 10~ 3 8.48 + 3x10^ 

1.000 8.72 ±3 x 10" a 8.73 ± 3 x 10" a 8.57 ± 2 x 10" a 



for Sl = 14. Fig. flW a) shows the dependence of the 

width er e = (\J (ej — (e.;)) 2 ^ for the onsite energy on 

the disorder strength Sd for Sl = 14. It can be seen 
from the graph that a increases linearly with the disorder 
strength, approximately equal to the speckle potential 
shift. Hence, the width is an appropriate measure of 
the disorder strength. The linear fit of this functional 
dependence gives <7 e = 0.95 x Sd- The distribution width 
of nearest neighbor hopping coefficients and Hubbard U 
are shown in Fig. [TST b) and Fig. fTBT c) respectively. In 
Fig. UGf d) , we show the disorder dependence of the mean 
value of Hubbard U. It can be seen from the graph that 
(u) is not sensitive to the disorder strength. 



V. CONCLUSION 

In this paper, we developed a method to construct low 
energy basis states and applied the method to calculate 
the Hubbard parameters in a disordered lattice created 
by an optical speckle field. The imaginary time projec- 
tion method introduced in this paper generates a type 
of Wannier-like localized basis that satisfies several con- 
ditions imposed by a reasonable coarse-grained, effective 
lattice Hamiltonian. 

Detailed many-body calculations using the deter- 
mined parameters with comparison to experiments are 
in progress. The method can be extended to include in- 
teractions. 
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